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In this paper, we study the Ising model with general spin S in presence of an external magnetic 
field by means of the equations of motion method and of the Green's function formalism. First, 
the model is shown to be isomorphic to a fermionic one constituted of 25 species of localized 
particles interacting via an intersite Coulomb interaction. Then, an exact solution is found, for any 
dimension, in terms of a finite, complete set of eigenoperators of the latter Hamiltonian and of the 
corresponding eigenenergies. This explicit knowledge makes possible writing exact expressions for 
the corresponding Green's function and correlation functions, which turn out to depend on a finite 
set of parameters to be self-consistently determined. Finally, we present an original procedure, based 
on algebraic constraints, to exactly fix these latter parameters in the case of dimension 1 and spin 
§. For this latter case and, just for comparison, for the cases of dimension 1 and spin 1 [F. Mancini, 
Eur. Phys. J. B 45, 497 (2005)] and spin 1 [F. Mancini, Eur. Phys. J. B 47, 527 (2005)], relevant 
properties such as magnetization (S) and square magnetic moment (S 2 ), susceptibility and specific 
heat are reported as functions of temperature and external magnetic field both for ferromagnetic and 
antiferromagnetic couplings. It is worth noticing the use we made of composite operators describing 
occupation transitions among the 3 species of localized particles and the related study of single, 
double and triple occupancy per site. 



I. INTRODUCTION 

In a previous article^, we have shown that a system 
constituted of q species of particles, satisfying Fermi 
statistics, subject to finite-range interactions and local- 
ized on the sites of a Bravais lattice, is exactly solvable 
in any dimension. By exactly solvable we mean that it is 
always possible finding a complete, finite set of eigenener- 
gies and eigenoperators of the Hamiltonian which closes 
the hierarchy of the equations of motion. Then, formal 
exact expressions for the corresponding Green's functions 
and correlation functions can be derived. The solution is 
only formal as these latter turn out to depend on a finite 
set ofparameters to be self-consistently determined. In 
Refs. IjlJ we have shown how is possible to exactly fix 
such parameters, by means of algebraic constraints, for 
the one-dimensional q = 1 and q = 2 cases, respectively. 

In Ref. U we have also shown that this system is iso- 
morphic to a spin- 1 Ising-like model in presence of an ex- 
ternal magnetic field. According to this, the exact knowl- 
edge of a complete set of eigenoperators and eigenener- 
gies of the system in any dimension acquires an evident 
relevance with respect to the hoary problem of solving 
the two-dimensional (in presence of an external magnetic 
field) and the three-dimensional Ising models. In partic- 
ular, the exact knowledge of the eigenenergies of the sys- 
tem can shed some light on the energy scales ruling the 
physical properties and the magnetic response of the sys- 
tem and can find an application as unbiased check for the 
approximate solutions present in the literature. More- 
over, this approach can eventually open a new route in 
the quest for an exact solution for these systems in higher 
dimensions as we have shown that it is always possible 
finding a formal exact expressions for the corresponding 
Green's functions and correlation functions. The solu- 



tion is only formal as these latter turn out to depend on 
a finite set of parameters to be self-consistently deter- 
mined. We have also shown how is possible to exactly fix 
such parameters, by means of algebraic constraints, for 
the one-dimensional case and we are now working on the 
possibility to use other algebraic constraints and topo- 
logical relations for higher dimensions. 

In this article, we apply this formulation to the one- 
dimensional q = 3 case. This model is isomorphic to 
the one-dimensional Ising model for spin S — | in pres- 
ence of an external magnetic field. This latter model 
can be exactly solved by means of the transfer matrix 
method which reduces the solution to an eigenvalue prob- 
lem of the fourth order (generally, we have an eigen- 
value problem of the 2S + 1-th order). Unfortunately, 
the analytic solutions of a fourth order equation are well 
known, but absolutely untractable and any result, along 
this way, could be obtained only by means of numeri- 
cal techniques, which surely do not facilitate a micro- 
scopic understanding of the properties of the system (see, 
for instance, Ref. In 1967, by means of a perturba- 
tion method, Suzuki, Tsujiyama and Katsura managed 
to reduce the order of the problem to the largest integer 
smaller than S + 1 in the case of zero external magnetic 
field and exactly computed the energy, the specific heat 
and susceptibility 5 . On year later, Obokata and Oguchi 6 
managed to apply the Bcthe approximation, which be- 
comes exact in one dimension, to the system and to re- 
cover the same exact results obtained by Suzuki. These 
latter were also recovered, few years later, by Silver and 
Frankel-:, who managed to formulate the problem as a dif- 
ference equation of order 25*+ 1 by an inductive technique 
and to reduce the order to the largest integer smaller than 
S + 1 in the case of zero external magnetic field. 

The spin- 1 Ising model in one dimension shows mag- 



2 



netic plateaus&SiiSiAi, i.e. topological quantization of the 
magnetization at the ground state of the system due to 
magnetic excitations, leading to the qualitatively same 
structures of the magnetization profiles of the Heisen- 
berg model. According to this, the study of this classical 
system can shed some light on the plateau mechanism 
and if it has purely quantum origin or can also depend 
on dimerization, frustration, single-ion anisotropy or pe- 
riodic field. These plateaus have been predicted not only 
in theoretical study but also have been observed in ex- 
perimental studies. For example, Narumi et alm^i^ ob- 
served the magnetic plateaus in the magnetization curve 
for both [Ni2(Medpt)2(n-ox)(H 2 0)2}(ClO i )22H 2 and 
[Ni 2 (Medpt) 2 {n - ox)(n - N 3 )}(Cl0 4 )o. 5 H20. Goto et 
alJA reported the existence of the magnetization plateau 
at 0.25 in spin-1 3,3,5,5-tetrakis (N-tert-butylaminxyl) 
biphenyl (BIP-TENO). In three dimensions and with 
additional terms in the Hamiltonian, the spin-| Ising 
model have been initially introduced to give a qualita- 
tive description of phase transition observed in the com- 
pound DyVOjp^ and also to describe ternary mixtures 16 . 
In particular, the phase diagram is not well known in 
contrast to the case S = 1, which represents a spe- 
cial case of the Blume-CapeliLiSiiiLSi and of the Blumc- 
Emery-Griffiths 21 models often used to study a vari- 
ety of interesting physical systems and, in particular, 
3 He — A He mixtures, fluid mixtures and critical phe- 
nomena. The mean field treatment^ predicts that the 
phase diagram differs for integer and half-odd-integer 
spins and does not present any multicritical point. While 
renormalization-group calculations 2 ^* 2 ^ and Monte Carlo 
simulations^ suggest the existence of a multicritical 
point, transfer matrix and conformal invariance studies^ 
show that there is no multicritical point in the phase di- 
agram. Unfortunately, Bcthe-PeierlsSZ^ and two-spin- 
clustcr approximations^ have not been performed in the 
low-temperature region. In qualitative agreement with 
the mean-field analysis, we can also find results from the 
self-consistent Ornstein-Zernike approximation^. 



II. THE MODEL 

We have analyzed a system constituted of q species of 
interacting particles obeying Fermi statistics, localized 
on the sites of a Bravais lattice and whose dynamics is 
ruled by the following grand-canonical Hamiltonian 



H = -n^2n a (i) + - ^2 Vab(hi)n a {i)n b (j) 



(2.1) 



ij, ab 



where i stands for the lattice vector R, and i = (i, t). 
n a {i) = <4(z)c a (z) ls the particle density operator of par- 
ticles of species a at the site i. c a (i) and c' a (i) are anni- 
hilation and creation operators, respectively, of particles 
of species a at the site i. They satisfy canonical anti- 
commutation relations 



{c a (i,f),4(j,i)} = SabSij 
{c a (i,t),c b (j,t)} = {cl(i,t),cl(j,t)} = 



(2.2) 



/i is the chemical potential and V a b(i,j) is the strength 
of the interaction between particles of species a and b at 
distance |i — j|. We have supposed that the particles are 
frozen on the lattice sites as their masses are very large 
and/or the interactions are so strong that the kinetic en- 
ergy is negligible. 

In this manuscript, we have focused on the case in 
which the interaction V b(i,j) does not depend on the 
particle species and is effective only between nearest- 
neighbor sites. That is, V a b(i,j) — 2dVa\j where d is the 
dimensionality of the system, ay is the projector on the 
nearest-neighbor sites and V is the bare strength of the 
interaction. For a hyper-cubic lattice of lattice constant 
a we have 



^ = ^E eik(R ^ Rj) «( k ) 

k 

1 d 

a(k) = - 2_j c °s(k n a) 



(2.3) 



In this manuscript, we present the exact solution of 
the model in presence of an external magnetic field. In 
the first section, we present the model. In the follow- 
ing section, we give the general solution in terms of the 
eigenoperators and of the eigenenergies. In the third sec- 
tion, we specialize the solution to the one-dimensional 
q = 3-S = | case. In the fourth section, we show how to 
close the self-consistent equations and compute all rele- 
vant correlation functions. In the fifth section, we present 
relevant properties such as magnetization (S) and square 
magnetic moment (S 2 ), susceptibility and specific heat as 
functions of the temperature and the external magnetic 
field both for ferromagnetic and antiferromagnetic cou- 
plings. Concluding remarks follow. 



where N is the number of sites. 

Hereafter, for a generic operator $(z), we will use the 
following notation 



$ Q (i,f)=5> sj $(j,t) 



(2.4) 



. It is also useful to introduce the vectorial notation 



c(i) 



fc 1 (i)\ 
c 2 (i) 



(2.5) 



4M 



C t(i)=( C t(<) 4(7) 

and to define the total particle density 
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Then, the Hamiltonian (|2.1|l becomes 



H = ~H^2 n(i) +dV^2 n(i)n a (i) 



Let us consider now the transformation 
n(i) = | + S(<) 

It is clear that 

n(i) = & S(i) = -q/2 
n(i) = 1 S(i) = 1-^/2 

n(i) = g-lo S{i) = q/2-l 
n(i) = q<=> S{i) = q/2 



(2.7) 



(2. 



(2.9) 



Under the transformation l]2.8[l . the Hamiltonian (|2.7() 
can be cast in the form 



H = -djJ2 s (. i ) s °'(. i )- h y^ s ( i ) + E o (2-10) 



where we defined 

J = -dV 
h = /! — g<iy 

£ = |(-M+f^)iV 



(2.11) 



The Hamiltonian (|2.10() is just the <i-dimensional spin-| 
Ising model with nearest-neighbor interaction in presence 
of an external magnetic field. In this manuscript, we have 
chosen to use the particle notation (|2.7|) . but the results 
that we have obtained are obviously valid for both the 
particle l|2.7|l and spin l|2.10|l systems after the transfor- 
mation l|2.8|) and the definitions (|2.11|) . 



III. COMPOSITE OPERATORS AND 
EQUATIONS OF MOTION 

In order to apply the equations of motion method and 
the Green's function formalism, we need to identify a 
suitable operatorial basis^i*^. It is immediate to ver- 
ify that the particle density operator n(i) has no time 
dependence 



i^n(i) = [n(i),H] = 



(3.1) 



According to this, the operator n(i), although it would 
have been a natural choice as component of the oper- 
atorial basis for such a system, is not suitable for this 
purpose. Let us introduce, instead, the following series 
of composite field operators 



p 



1,2,. 



(3.2) 



whose first element is just c(i). These fields satisfy the 
following hierarchy of equations of motion 



-Hi/j p (i) + 2dVip p+ i{i) 



(3.3) 



Now, because of the anti-commutation relations (|2.2|l . 
it can be shown that the operators [n a (i)] p satisfy the 
following relation 



2qd 



(3.4) 



where the coefficients i„ are rational numbers that can 
be easily determined after the algebra and the actual 
structure of the lattice (see Appendix A for d = 1 and q = 
3). Then, for p — 2qd + 1, the hierarchy of equations of 
motion l|3.3|l closes as the additionally generated operator 
ip2qd+2(i) = c(i)[n a (i)] 2qd+1 can be rewritten in terms of 
the first 2qd + 1 elements of Q3.2H through the relation 
(El- 
According to this, the n-component composite field op- 
erator ?/>(i), defined as 



( 



c(z) 
c(i)n a (i) 



\ c(i)[n a (i)] r 



\ 
I 



(3.5) 



where n — 2qd+l, is an eigenoperator of the Hamiltonian 
(|2.7|) as it satisfies the equation of motion 



iJ^(*) = [lK*)»#]=^(*) 



(3.6) 



where the n x n energy matrix e has the following ex- 
pression 







2dV 
-A* 





2dV 



2dV 



\ 



V 2dVA[ n) 2dVA ( ™ ] 



- f i + 2dvA i l ; i 2 1 j 

(3.7) 



The eigenvalues E n of the energy matrix e read as 

E m = -fj, + (m - 1)V m = l,2,...,n (3.8) 

After (|3.6|l . we can claim to have formally, but exactly, 
solved both Hamiltonians (|2.7|l and l|2.10|l as we have de- 
termined for them a complete set of eigenoperators and 
eigenvalues in any dimension d. The solution is only for- 
mal as we have still to compute all correlation functions. 
Let us define the thermal retarded Green's function 

G R (i,j) = {Rim^Hm = o{u-t 3 )({^{i),^{])}) 

(3.9) 
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where (• ■ • ) denotes the quantum-statistical average over 
the grand canonical ensemble. By introducing the 
Fourier transform 



G R (i,j) 



^ ff dLsdke fk -&*-' B *)- iu < tt - t AG R (k,w) 



(3.10) 



and by means of the equation .(j|> , we obtain the equa- 
tion 



(3.11) 



where J(k) is the Fourier transform of the normalization 
matrix 



7(i,j) = ({V(i,t)> + (j,t)}) 
Equation 13.111) gives 

^ w - £m + i<5 



(3.12) 



(3.13) 



The spectral density matrices a^™^ (k) can be calculated 
by means of the formula™ 



(3.14) 



where f2 is the nxn matrix whose columns are the eigen- 
vectors of the matrix e. The correlation functions 

C(i,j) = (ip(i)^(j)) (3.15) 

can be immediately calculated after (|3.13|) and read as 



where 



m— 1 



r,„ = I + iaiil) ( 



(3.16) 



(3.17) 



T is the temperature. By similar techniques we 
can calculate multi-point correlation functions as 
C(i,j;h,l 2 ,...,l s ) = mW(j)n{h)n{l 2 )---n{l s )). For 
an example, see Ref. 

Formally, (|3.13|l and (|3.16|) constitute the exact solu- 
tion of the problem under analysis. The solution is only 
formal as the complete knowledge of the Green's func- 
tion and of the correlation functions is not fully achieved 
owing to unknown static correlation functions appearing 
in the normalization matrix J(k) because of the complex 
non-canonical algebra satisfied by the field ip(i). The 
unknown correlators are expectation values of operators 
not belonging to the chosen basis tp(i) and should be 
self-consistently calculated in order to definitively con- 
clude the analysis and get the complete exact solution. 



Unfortunately, the derivation of a set of self-consistent 
equations capable to determine the unknown correlators 
is not an easy task at all. Moreover, it depends on the 
particular model under analysis and, in particular, on its 
dimension d. In Ref. |l|, we have developed a technique to 
obtain a complete exact solution in the one-dimensional 
case (d = 1) by determining the unknown correlators by 
means of algebraic constraints. Within this framework, 
we have studied the 2-statc 2 (q = 1) and the 3-state 3 - 
(q = 2) models. In this article, we consider the one- 
dimensional 4-state (q = 3) model. 



IV. THE ONE-DIMENSIONAL 4-STATE MODEL 

We now apply the general formulation given in the pre- 
vious Section to an infinite homogeneous chain (d = 1) 
with q = 3. The derivation of the relation l|3.4|l for this 
specific case and the explicit expressions of the coeffi- 
cients 4m are given in Appendix A. Accordingly, the 
last row of the energy matrix e has the following non- 
zero entries 



£72 



£75 



45 441 



£74 = -203U 



735 



V £76 



175 



(4.1) 



V e 77 = -n + 2lV 



The matrix f2 has the expression 



= 



/I 







Vo 



2 6 

2 5 
2 4 
2 3 
2 2 
2 
1 



(2/3) 6 
(2/3) 5 
(2/3) 4 
(2/3) 3 
(2/3) 2 
(2/3) 
1 



(1/2) 6 
(1/2) 5 
(1/2) 4 
(1/2) 3 
(1/2) 2 
(1/2) 
1 



(2/5) 6 
(2/5) 5 
(2/5) 4 
(2/5) 3 
(2/5) 2 
(2/5) 
1 



(V3) 6 \ 
(1/3) 5 
(1/3) 4 
(1/3) 3 
(1/3) 2 
(1/3) 
1 / 



(4.2) 



and the normalization matrix / reads as 



(hi 


hi 


h 3 


hi 


Jl5 


h6 


Jl7 


\ 


hi 


h 3 


In 


Jib 


Im 


In 


J27 


ha 


hi 


hs 


Jie 


In 


I 27 


J37 




hi 


h 5 


he 


In 


I 27 


I 37 


J47 




h 5 


he 


In 


In 


I 37 


J47 


J57 




he 


In 


In 


J37 


J47 


J57 


J 6 7 




\Il7 


In 


In 


J47 


J57 


J 6 7 


J77 


/ 



(4.3) 



where 

h, P = Kb- 1 ) 



Ip,7 = £ h, m+ iA%+V 



= ([n a (i)f> 



(4.4) 
(4.5) 
(4.6) 



Then, the spectral density matrices <r^ m ' can be easily 
computed by means of Eq. (|3.14|) and the correlation 
functions by means of Eq. I|3.17|l . 
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It is worth noticing that, according to the struc- 
ture of the normalization matrix I [ 14.3( 1 and (|4.5I) ]. 
which is dictated by the relation ((3.4(1 . there exist only 

seven independent spectral density matrices cr%^ and 
seven independent two-point correlation functions C\, v = 
(c(i)c^(i)[n Q (z)] p_1 ). All others can be obtained as com- 
binations of these latter according to 1(3.4(1 and 1(4.5(1 . 
As regards the spectral density matrices, we have 



correlators are unknown as the operators [^ Q («)] p do 
not belong to the chosen basis 13.511 . 

To the purpose of definitively computing the correla- 
tion functions Ci,fc, let us introduce the following opera- 
tors 



r (m) 



(m) 



(4.7) 



£,(*) = [l-n(t)+-D(i)]c„(i) 
T) a {i) = [n(i)-2D{i)]c a (i) 
C a {i) = D(i)c a (i) 



(5.1) 



where S m are functions of the elements Ji, p with p = 
1, . . . , 7 and 1^™) are numerical matrices. In particular, 
we have the expressions 



where the double D(i) and the triple T(i) occupancy op- 
erators at the site i read as 



Si = Ji,i 

+ i-(-441/i, 2 + 8127i,3 



7357i, 4 + 350/i,5 - 84Ji, 6 + 8/i, 



a ,6=1 
a<b 

3 

T(t)=n».w 



(5.2) 



£ 2 = — (90/i,2 - 261/i,3 + 290/1,4 - 155/1,5 + 40/ 1>6 - 4/i, 7 ) o=i 

£ 3 = i (-90/i,2 +351/i,3 -461/i,4 + 274/i,s - 76/i, 6 + 8/1,7) The operators £„(£), rj a (i) and ( a (i) rule the transitions 
" among states with different particle number at the site i: 

£ 4 = I (30/i,2 - 127/i,3 + 186/i,4 - 121/i,5 + 36/i,6 - 4/1,7) 1, 1 2, 2 3, respectively. Their complex algebra 

is presented in detail in Appendix B. In particular, they 



£5 
£ 6 



(-45/i,2 + 198/i,3 - 307/i,4 + 214/i,5 - 68/1,6 + 8/1,7) satisfy the following relevant relations 



(18/l,2 - 81/l,3 + 130/l,4 - 95/l,5 + 32/i,6 - 4/l,7) 



£7 = —(-30/i,2 + 137/1.3 - 225/i,4 + 170/i,5 - 60/i,6 + 8/1,7) 
yu 



(4.8) 



Ca(i) = &,(*) +Va(i) +Ca(«) 

and 



(5.3) 



and 



p(l) 

1 l,m 


= (1 0) 








r (2) 
1 l,m 


= (1 2" 1 2" 2 2" 3 2~ 4 2" 5 2 


- 6 ) 






p(3) 
1 l,m 


= (1111111) 








r (4) 
1 l,m 


= (1 (2/3)" 1 (2/3)" 2 (2/3)" 


3 (2/3) " 4 (2/3) " 5 


(2/3)" 6 ) 


for any 


p(5) 
1 l,m 


= (1 (1/2)- 1 (1/2)- 2 (1/2)" 


3 (1/2)- 4 (1/2)- 5 


(1/2)" 6 ) 


After 
C) 


r (6) 

1 l.m 


= (1 (2/5)- 1 (2/5)- 2 (2/5)- 


3 (2/5)- 4 (2/5)- 5 


(2/5)- 6 ) 


p(7) 
1 l.m 


= (1 (1/3)" 1 (1/3)" 2 (1/3)- 


3 (1/3)- 4 (l/3)~ 5 


(1/3)" 6 ) 





nUi)D(i) = 



eI(i)T(t) = 
»?t(i)r(i) = 

Ci(W) = o 

(5.4) 



After (|5.4|) . it is possible to show that (see Appendix 



-f3H 



-0H oi 



(4.9) 

As regards the two-point correlation function, we have 
1 I 

C(i,j) = S i} - T m a^e- iE ^-^ (4.10) 

m— 1 

V. SELF-CONSISTENT EQUATIONS 

The correlation functions depend, through the 

spectral density matrices that are functions of the en- 
tries of the normalization matrix 1(3.14(1 . on the correla- 
tors = ([n a (i)] p ) with p = 1, . . . , 6. At this stage, the 



C\i)e 



-f3H 



C+(i)e 



^(i)^2(2f m +g m )[n a (i)] m e 



-/3H 0i 



m-/3H 0i 



where Hoi is a reduced part of the Hamiltonian 

Hoi = H — Hji 
Hn = 2Vn(i)n a (i) 



(5.5) 



(5.6) 



We have used for £(z), rj(i) and ((i) the same vectorial 
notation that has been used for c(i) (see l(2.5(l ). The 
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coefficients f m and g m are defined in Appendix C. We 
will exploit the relations (|5.5|l in order to compute the 
correlation functions 



C^ = (v(^(i)[n a ^} k - 1 ) 
C C 1 i = (t(i)<; i (i)[n a (i)] k - 1 ) 



k = 1, 



(5.7) 



that, after the relations l|5.3[l . sum up to C± t k = Cf 
Ox\ + W k- By means of (|5.5J) . we have 



ii 



« _ (e(Qg f (0[n a (Q] fc - 1 )oi 

k (e-^«) w 



„„ _ (7 ? (») ?? t(z)K(z)] fc - 1 ) 0t 



(e-^ ff «>(H 

(r 7 (z)7?t(z)[n"( ? )] fc - 1 +"') 0i 



m— 1 



<e-^>o; 



(5.8) 



where for a generic operator O the notation (O)oi denotes 
the thermal average with respect to Hqi 



(O)oi 



Tr{Oe- fJHo >} 
Tr{e-P H °>} 



(5.9) 



Now, we observe that H$i describes a system where 
the site i is not connected to any other site of the chain. 
Then, at the site i the local operators enjoy a free dy- 
namics 



W),Hoi] 






[r](i),H 0i \ 




(5.10) 




= -K(<) 





By means of these equations of motion and by making 
use of the relations JBSJl, it is possible to derive 



(f»r(<)[« a (or l ) w = (1+egM)3 <[" a (0] fc - 1 ) M 

Op/3 A 1 

(5.11) 



and rewrite l|5.8[l as 

i (K( i )] fe - 1 ) 0i 



r<nn 



(l + e^)3 ( e -/Jff«) w 
2e^ (K(z)] fe - 1 } 0l 



(l + e^)3 ( e -/3H») 0i 

2e^ ^ . ([n a (i)p+ fc -% 



(1 + 6/3^)3 ^ 
x ' m— 1 



([n a (t)] fc_1 )w 



(l + e^)3 ( e -/Jff«) w 
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(1 + e^)3 ^ 

7 m— 1 



(2/ m +3m)' 



[n«(*)] 



m+fc— 1 



(5.12) 



Therefore, we are left with the problem of computing the 
functions {[n a {i)] p ) ov 

Now, we observe that describes a system where the 
original infinite chain is split into two disconnected infi- 
nite sub-chains (the infinite chains to the left and to the 
right of the site i) plus the site i. Then, in the H^- 
representation, correlation functions which relate sites 
belonging to different sub-chains and/or the site i, can 
be decoupled: 



(a(j)b(l)) 0i = (a(j)) 0i (b(l)) 0i 



(5.13) 



when, for instance, j < i < 1. a(j) and b(l) are any 
functions of n(j) and n{l), respectively. Let us recall 
Eqs. ijA7f> and (|A8I) . By using the property l|5.13f) . we 

have 



(5.14) 



Zo(i))oi 


= 2X X 


Zi(i)) 0i 


= 2X 2 + X\ 


Z 2 (j)) 0l 


= 2X3 + 2X1X2 


Z 3 (i))oi 


— 2X1X3 + X2 




= 2X2X3 


Zs{i))oi 





with 




Xi = 




x 2 = 


(£>«(*)> 


x 3 = 


(T Q (i)) ( 



(5.15) 



Therefore, we have 
1 



(K(i)n« = ^ 2X 1+ b[ p \2X 2 + Xl) 



-2b^(X 3 + 
-2^^X2X3 



X 1 X 2 ) + b 3 p) (2X 1 X 3 + Xl) 



(5.16) 



b^xl 



The problem of computing all two-point correlation func- 
tions is thus reduced to the problem of computing just 
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three parameters: Xx, X 2 and X 3 . If we suppose the 
system to be homogenous, we can use the following three 
self-consistent equations in order to determine the un- 
known parameters 



(n(<)> = (n a (i)) => (n(*)e-^ l )«H - (n a (i)e-^) 0i 
(D(i)) = (D a (i)) (D(i)e-^) 0i = {D a {i)e~^) Qi 
(T(i)> = (T a (i)) =s> (T(*)e-^) w = (T a (i)e-^) 0i 

(5.17) 

By means of Eqs. I|C2|I and (|B9|I . we have 

6 

( n (i)e-^«) w = Si + E + 2B2 ) 

m— 1 

+.g m (25 2 + 3B 3 ) + ^ m 3B 3 ] <[n Q (i)] m > 0i 

8 (5-18) 

(7i a (i)e-^»> . = X x + [f m Bx+g m B 2 



+h m B 3 ]([n a (i)} m+1 ) 



{D(i)e-P Hli ) oi = B 2 +J2 [(25a + 3S 3 )/„ 

m— 1 

+ (5 2 + 6B 3 )9m + 3B 3 h m ] {[n a (i)] m ) oi 

6 

(^ a (z)e- /3H ^)o J = X 2 + [Bifm + B 2 g ri 

m—1 

+B 3 h m ] (D a (i)[n a (i)] m )oi 

6 

{T(i)e-^) 0l =B 3 ( 3 /m + 35m 

m—1 

+h m )\ {[n a (i)] m ) oi 

6 

{T a (i)e- pH ^) ai = X 3 +Y, {Bxf m + B 29n 

m—1 

+B 3 h m ) {T a {i)[n a (i)] m ) ai 
where 

Bx = (n(i))oi = 
S 2 = (£>(i))oi = 



B 3 = (T(i))« = 



e /3M + 1 

3e 2 ^ 

( e /3M+ 1)2 

(e"" + l) 3 



(5.19) 



(5.20) 



(5.21) 



We need to calculate the averages (D a (i)[n a (i)] m )oi 
and (T a (i)[n a (i)} m ) 0l . By using (£3 



m=0 
5 



(5.22) 



<T«(i)[n Q (i)f) 0i = &«<T«(i)Z m (i)) 



By recalling the definitions (|A8|) and using the property 

(£>«(*) [n^Hex = ^ [2X 2 + 3X 3 + X!X 2 
+b[ p) {X 2 + 6X 3 + Xl + 2X X X 2 + 3X X X 3 ) 
+b { 2 p) (3X 3 + 2Xl + XiX 2 + 6XxX 3 + 4X 2 X 3 ) 
+b {p) (3XiX 3 + Xl + 8X 2 X 3 + 3Xl) 



+2b {p) (2X 2 X 3 + 3X 3 2 ) + 3b {p) X 2 3 



(5.23) 



{T a {i)[n a (i)Y) 0i = -[ZX 3 +XxX 3 
+b[ p) {3X 3 + 3X 1 X 3 + X 2 X 3 ) 

+b (p) (X 3 + Xl + 3X2X3 + 3X X X 3 ) (5-24) 
+6 3 p) (3X 2 X 3 + X X X 3 + 3X 3 2 ) 
+6i p) (3X 3 2 + X 2 X 3 )+4 p) X 3 2 



Summarizing, the three parameters Xx, X 2 , X 3 
are determined by the coupled self-consistent equa- 
tions ijETTgjl - 1(530)1 . where the averages ([n a (i)] m ) oi , 
(D a (i)[n a (i)] m ) Q i and (T a (i )[n a ( i)] m ) 0i are computed by 
means of (|5.16|l . (|5.23() and (|5.24(l . respectively. Once we 
know the three parameters, we can calculate the correla- 
tion functions and all the properties. For example, 



"" - 2Bl ~ 4Ba + 6Ba [([n a (f)] k - 1 ) [ 



°l,fc - 



3(e-^)c 



3(e-^)c 



£(2/ m +5 m )(K(i)] m+fc - 1 > J 



Cx,k 


= eg, + c\\ 




(n) 


= 3 — 3Ci ; /c 




(D) 


= 3 - 3(C5 + 


|cf5 1 + 2C^) 


(T) 


= 1-3(±C« 


+ 5^1+^1) 









(5.25) 



+ 2_^\ t >XJm + 2 g m +U 3 tl m ) — - _n H , — 
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(n(i)[n a (i)] p ) = B x 

6 



[n a {i)] p ) 



Oi 



e HO 

[(Bi + 2B 2 )f m + {2B 2 + 3B 3 )g„ 

m—1 



(5.26) 



+3B 3 h. 



-PHl) 



The coefficients f m , g m and h m are defined in Appendix 
C. The average can be computed by means of 

(|C2() and has the expression 

6 

m—1 

(5.27) 

VI. RESULTS 
A. Ferromagnetic coupling J = 1 

In Fig. ^ we report the magnetization per site m as 
a function of the magnetic field h at J = 1, the double 
D and the triple T occupancies per site as functions of 
the chemical potential /i at V = — 1 , the square magnetic 
moment (S 2 ) per site as function of the magnetic field h 
at J = 1 and the magnetic susceptibility x as a function 
of the magnetic field h at J = 1 and T = (only in mid- 
dle bottom panel), 0.5 (except for middle bottom panel), 
1 and 2. At the lower temperature (T = 0.5), going from 
negative to positive values of the magnetic field, the mag- 
netization jumps from — | to | and the susceptibility di- 
verges exactly at zero field. This is the signature of a 
quite sharp transition between spin configurations with 
maximum possible magnetization and direction dictated 
by the value of the external magnetic field. This immedi- 
ately reflects on the values of the double and triple occu- 
pancies that quite sharply jumps from zero to their max- 
imum possible values: 3 and 1, respectively. The tem- 
perature just makes the transition smoother and forces 
the system to pass through excited states with a lower 
value of the absolute magnetization as it results clear by 
looking at the values of the square magnetic moment. 

In Fig. |2 we report the inverse magnetic susceptibil- 
ity x _1 at h — and the magnetic susceptibility \ a t 
ft, = 0.1 as functions of the temperature T at J = 1 and 
spin 1/2, 1 and 3/2. The sharp transition at zero tem- 
perature in zero external magnetic field is here clearly 
signalled by the divergence of the magnetic susceptibility 
for any value of the spin. At finite external magnetic field 
h = 0.1, the susceptibility presents a maximum at a fi- 
nite temperature whose value depends on the value of the 
spin. Both at very low and very high temperatures, the 
magnetization stays constant at its maximum value or at 
zero, respectively, almost independently on the value of 
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FIG. 1: (top) Magnetization per site m as a function of the 
magnetic field h at J — 1. (middle top) Double D and triple 
T occupancies per site as functions of the chemical potential fi 
at V = — 1. (middle bottom) Square magnetic moment (S 2 ) 
per site as function of the magnetic field h at J = 1. (bottom) 
Magnetic susceptibility \ as a function of the magnetic field 
h at J = 1. In all panels, T — (only in middle bottom), 0.5 
(except for middle bottom), 1 and 2. 
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FIG. 2: (top) Inverse magnetic susceptibility x _1 at ft = and 
(bottom) magnetic susceptibility x at ft = 0.1 as functions of 
the temperature T at J = 1 and spin 1/2, 1 and 3/2. 



3.5 




FIG. 3: Specific heat C as a function of the temperature T 
at J = 1 and h = 0.1, 0.5, 1 and 1.5. 

the external magnetic field. According to this, in these 
two regimes the susceptibility is zero or practically zero. 
At intermediate temperatures, the value of the magne- 
tization definitely depends on the value of the external 
magnetic field. This explains the bell-shape of the sus- 
ceptibility curve and the presence of a maximum. 

In Fig. EI we report the specific heat C as a function of 
the temperature T at J = 1 and h — 0.1, 0.5, 1 and 1.5. 
The exponential behavior at low temperatures clearly 



signals the presence of at least one gap in the excita- 
tion spectrum between the fully polarized state, which is 
the ground state for any non-zero value of the external 
magnetic field, and the first excited state to which corre- 
sponds a lower value of the magnetization. This behavior 
also explains why we have been able to see a sharp tran- 
sition in the magnetization as a function of the external 
magnetic field at temperatures as high as T — 0.5. The 
temperature Tm, at which is clearly observable a maxi- 
mum in all curves, strongly depends on the value of the 
external magnetic field h. In particular, Tm increases on 
increasing h: the gap between the ground state and the 
first excited state also increase on increasing h. We can 
speculate that Tm ~ h(h + 2J), that is half of the gap 
existing between the fully polarized states with spin | 
and |. 



B. Antiferromagnetic coupling J = — 1 

In Fig. 0] we report the magnetization per site m as a 
function of the magnetic field h at J = — 1, the square 
magnetic moment (S 2 ) per site as function of the mag- 
netic field h at J = — 1 and the double D and the triple 
T occupancies per site as functions of the chemical po- 
tential \i at V = 1 and T = 0, 1 and 2. At zero tem- 
perature, going from negative to positive values of the 
magnetic field, the magnetization first jumps from — § 
to at h = 3J and then from to | at h = — 3J. At 
these values of the external magnetic field, the state of 
lower energy changes from the fully polarized one with 
spin — | (E = 0) to the antiferromagnetically aligned one 
with spins | and — I (E = — J) and from this latter 
to the the fully polarized one with spin ~ {E = — 3h), re- 
spectively. The double and the triple occupancies behave 
accordingly. The temperature just makes the transition 
smoother as it results clear by looking at the values of 
the square magnetic moment. These latter also show 
the typical antiferromagnetic behavior that will be easily 
identified in the features of the susceptibility discussed 
below. 

In Fig. [SI we report the magnetic susceptibility x a t 
J = —1 as a function of the external magnetic field h for 
T = 0.8, 1, 1.5, 2, 3 and the temperature T for h = 
and spin 1/2, 1 and 3/2. As we could expect from the 
magnetization curves, the susceptibility diverges, at zero 
temperature, at the critical values discussed above and is 
null for all other values. On increasing the temperature, 
the susceptibility becomes finite for any value of the ex- 
ternal magnetic field, although not negligible only in an 
increasing, but finite, range of values. According to this, 
as the integral of the curves over the whole axis should 
give a constant value, 3, that is, the larger jump in the 
magnetization, the heights of the two peaks steadily de- 
crease on increasing the temperature. As a function of 
the temperature, the susceptibility has the typical ther- 
mal activated behavior of the parallel antiferromagnetic 
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FIG. 4: (top) Magnetization per site m as a function of the 
magnetic field h at J = — 1. (middle) Double 7? and triple T 
occupancies per site as functions of the chemical potential fi 
at V = 1. In both panels, T = 0, 1 and 2. (bottom) Square 
magnetic moment (iS 2 ) per site as function of the magnetic 
field ft at J = — 1. 



susceptibility for any value of the spin. According to 
this, we can interpret the position of the maximum as 
the temperature at which the system loses any memory 
of its antifcrromagnetic ground state at zero temperature. 

In Fig. HO we report the specific heat C as a function 
of the temperature T at J = — 1 and h = — 5 => —1.5. 
The curves seem to show the typical pattern caused by 
the interplay between at least two gaps in the excitation 
spectrum: h — 3J and h — 2 J. These would correspond 
to transitions between the fully polarized state with spin 
— | and the antiferromagnetically aligned state with spins 





— T=0.8 




■■- T=lj0 




--1=15 




--1=2.0 




--T=3.0 




FIG. 5: Magnetic susceptibility x at J = —1 as a function of 
(top) external magnetic field h for T = 0.8, 1, 1.5, 2, 3 and 
(bottom) temperature T for h = and spin 1/2, 1 and 3/2. 



and — | 



and the fully polarized state with spin — b, 
respectively. The vanishing value of the first gap, in the 
range of values of the external magnetic field presented 
in the figures, could explain the much larger sensitivity of 
the magnetization to temperature with respect to what 
happens in the ferromagnetic state. It is worth noticing 
the appearance of crossing points^ for 3J < h < 2.5J 
and 2.5J < h < 1.5J. 



VII. CONCLUSIONS 

We have studied the Ising model with general spin S 
in presence of an external magnetic field by means of the 
equations of motion method and of the Green's function 
formalism. First, the model has been shown to be isomor- 
phic to a fermionic one constituted of 2S species of local- 
ized particles interacting via an intersite Coulomb inter- 
action. Then, an exact solution has been found, for any 
dimension, in terms of a finite, complete set of eigenoper- 
ators of the latter Hamiltonian and of the corresponding 
eigenenergies. This explicit knowledge has made possible 
writing exact expressions for the corresponding Green's 
function and correlation functions, which turn out to de- 
pend on a finite set of parameters to be self-consistently 
determined. Finally, we have presented an original proce- 
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— h=-2.9 




— h=-2.8 




--h=-2.1 




--h=-2.6 




- - h=-2.5 



c 1.0 




dure, based on algebraic constraints, to exactly fix these 
latter parameters in the case of dimension 1 and spin |. 
For this latter case and, just for comparison, for the cases 
of dimension 1 and spin \ and spin 1 , relevant properties 
such as magnetization (S) and square magnetic moments 
(S 2 ), susceptibility and specific heat are reported as func- 
tions of temperature and external magnetic field both for 
ferromagnetic and antiferromagnetic couplings. Ground 
state properties and relevant transitions and gaps have 
been studied. Crossing points in the specific heat have 
been identified. 



FIG. 6: Specific heat C as a function of the temperature T 
at J = — 1 and h = — 5 =>• —1.5. 



APPENDIX A: ALGEBRAIC RELATIONS 
1. Formula for n p (i) 

We want to find a recurrence relation for the operator n p (i), with p > 1. At first we note that 

n v(i) = K + n 2 + n 3 r = ELo ( l) ELo ( ? ) n?- n n%- l n l 3 = [n? + n\ + n\\ 

f \ V J V J ( \ ( \ ( A1 ) 



Because of the property n p = rife, IjAlfl takes the form 

n=l v 7 n=2 v 7 ;=1 v 7 

where and T(i) are the double and triple occupancy operators as defined in l|5.2J) . 
The sums in <|A2ll can be analytically performed 




= 3(l-2P + 3 p - 1 ) 



and we have the algebraic relation 
n p (i) = n(i) + b[ p) D(i) + b { i ] T{i) 

2. Formula for \n a (i)f 

By recalling that 
K[i)]=\[n{i + a)+n{i-a)] 
we have for p = 2, 3 

K«F - w E^= ( P m ) ™ p - m (* + «K l (* - a) 

= ±n?(i + a) + - a) + & E£~=i ( ^ ) " p_m (i + «)» m (* - «) 

Because of the algebraic relation 1A4|I we obtain 

where the operators Z m (i) are defined as 





= 2n a {i) 




Z x {%) 


= 2D a (i) + n(i 


+■ a)n(i — a) 


Z 2 (i) 


= 2T a (i) + n(i- 


f a)D(i — a) + n(i — a)D(i + a) 




— n(i + a)T{i — 


a) + D(i + a)D(i - a) + n(i - a)T(i + a) 


Z 4 (i) 


= D(i + a)T{i- 


-a)+D(i-a)T(i + a) 




= T(i + a)T(i- 


a) 



and the new coefficients bm have the expressions 
= 1 

&3 W = ££i ( J, ) bt" = 4(-l + 2 2 ^ 2 + 3 ■ 2^ - 3P) 
bi P) = Emil ( jjj ) b[ p - m) b { 2 m) = 5(1 - 2 2 *> - 2? +1 + 2 • 3P - 5^) 
4 P) = Emil ( J, ) 4 ? " m) &2 m) = 6(-l + 5 • 2P- 1 + 5 • 2 2 ?- 1 - 10 • 3P" 1 -5P + 6?- 1 ) 
By solving (|A7|) with respect to the quantities Z m (i) we obtain the recurrence relation 
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where the coefficients i„ are rational numbers defined as 



A? 




6 (p) 4- 2 

u l '3 


a 2 2 






2^ p) + 


6 °3 




~~ 2J> L 3 u 2 


26^ + 


3°4 




— 2p i 3 "3 


§4 P) + 


9 °5 J 




- J_[J-/>) , 

2P L 15 "4 


- W\ 




A (p) 


1 4 yip) 

2p 45 5 







1^4 



(P) 



5 dp) , 137 tW 



3 "4 
Pf 



+ 



90 °5 



(All) 



We note that 



.4 



(V) 



E 6 

A$=6 mp (p<6) 



(A12) 



In table 1 we give some values of the coefficients (p > 



p 


aw 


4 p) 


M 


4 P) 


A (P) 


A(P) 
A 5 


7 


45 


441 


203 


735 


175 


21 


4 


8 


2 


8 


4 


2 


8 


945 


9081 


8085 


13811 


735 


133 




8 


lfi 


8 


lfi 


2 


2 


9 


5985 


56763 


98915 


81585 


32739 


1323 




lfi 


lfi 


lfi 


lfi 


4 


10 


59535 


559503 


480375 


774575 


37485 


22827 


Ifi 


32 


lfi 


32 


4 


lfi 



APPENDIX B: ALGEBRA OF THE PROJECTION OPERATORS 

On the basis of the canonical anti-commutation relations 12.2fl it is straightforward to derive the algebra satisfied 
by the operators £ a (*), Va(i), Ca(*) defined in (|4.2|) . Their anti-commutation relations are 



(l-n 2 (i))(l-n 3 (i)) -ci(*)4(*)(l-n 3 (<)) -ci(*)4(*)(l - n 2 (i)) 

{£(*U + 0')} = *U | -09(t)(4(i)(l-n3(i)) (*))(! -*»3(i)) -oa(i)4(*)(l ~ 

-c 3 (i)4(i)(l-n 2 (i)) -c 3 (04(i)(l-m(i)) (1 -m(*))(l -na(i)) 



(Bl) 



( n 2 (i) + n 3 (i) - 2n 2 (i)n 3 (i) ci(i)4(i)(l - 2n 3 (i)) ci(i)4(«)(l - 2n 2 (i)) 

c 2 {i)c\{i){\ - 2n 3 (i)) m(i) + n 3 («) - 2n 1 (i)n 3 (i) c 2 (i)4(*)( 1 ~ 2ni(i)) 
\ c 3 (i)cj(i)(l - 2n 2 (i)) c 3 (i)4(«)(l - 2m(i)) m(i) + n 2 (i) - 2rn(i)n 2 (i) 



(B2) 



/ n 2 (i)n 3 (i) d(i)cl(i)n 3 (i) ci(i)4(i)n 2 (i) \ 
c 2 (i)4(i)n 3 (i) n 3 (i)ni(i) c 2 (i)4(«) n i(«) 



\ c 3 (i)c\(i)n 2 (i) c 3 {i)c\{i)ni{i) ni(i)n 2 (i) J 



{mxHj)} = h 

It easy to verify that 

{^U f (i)} + {#WC?)} + {C«,C + (i)} = *u 

Other relevant algebraic properties are 

ELi 6»(i)£(0 = 3(1 - + !>(») - T(i)) 
EL 1 ^a(<)=«(<)-2I>(*) + 3T(i) 



(B3) 




(B4) 



(B5) 



Etl f«(0^(0 = M*) - + 6T « 

ELi^W%(i) = 2 J D(i)-6T(i) 



(B6) 
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E!UC»(0Cj(i) = tf(i)-3T(i) 
ELiClWC«(<) = 3T(<) 



(B7) 



By means of these relations we can express the particle density operator, the double and triple occupancy operators 
as 

nii) = 3 - Ea=ite-0£l(i) + + C«WClW) 

= 3 - ELi(f-W^(i) + + 2Ca(*)Cl(i)) ( B8 ) 

= i - ELi(|€«(i)^(0 + ^a(i)»/i(i) + Ca(0Cl(0) 

Also, the following relations hold 

n 2 {i) = n(i) + 2D(i) D 2 (i) = D(i) + 6T(i) 

n(i)D(i) = 2D{i) + 3T(<) D(i)T(i) = 3T(i) (B9) 

n(/)T(i) = 3T(i) T 2 (z) = T(i) 



APPENDIX C: CALCULATION OF e~ pHl ' 



From the definition ifjj = 2Vn(i)n a (z)we have 

e - ^" = 1 + V -(-l) p (2f3V) p n p (i)[n a (i)] p (CI) 
p=i 

Recalling (|A4|) and (|A10|) we can write 



e-^« = 1 + E£=i ^-WWWO + &?°I>(i) + «#°T(i)] EL=i ^ [n a (i)P (C2) 
= 1 + nii) El =1 /m[n a (i)] ro + £l=i 5m K(*)P + r(i) EL=i M™ q (*)P 

where 

ffm=Er=iF(- i ) p (w p & i p) ^ ) (C3) 

By using the explicit expression (|A11|I of the coefficients , the infinite sums in l|C3(l can be analytically performed. 
Straightforward calculations show that the parameters f m , g m and h m are linear combinations of e~ n @ v with n ranging 
from to 18. ^ 
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